K <- function(u) {  
  if (any(u < -1)) 0  
  else {
    u[u > 1] <- 1
    (-1/4)^(length(u)) * prod((u+1)^2 * (u-2))  
  }
}

H <- function(data, x, h) {
  Kn <- array(0, c(dim(data)[1]))
  for (i in 1:dim(data)[1]) {
    u <- (x - data[i,]) / h  
    Kn[i] <- K(u)
  }
  (1/dim(d)[1]) * sum(Kn)
}

H2 <- function(data, x, l, h) {
  n <- dim(data)[1]
  d <- dim(data)[2]
  Kn <- array(0, c(n))
  for (i in 1:n) {
    u.upper <- (x + l - data[i,]) / (h * data[i,])
    u.lower <- (x - data[i,]) / (h * data[i,])      
    if (any(u.upper < -1, u.lower > 1)) Kn[i] <- 0
    else {      
      u.upper[u.upper > 1] <- 1
      u.lower[u.lower < -1] <- -1            
      Kn[i] <- prod(u.upper - 1/3 * u.upper^3 - u.lower + 1/3 * u.lower^3)    
    }
  }
  (1/n) * (3/4)^d * sum(Kn)  
}

H3 <- function(data, x, l, h) {
  n <- dim(data)[1]
  d <- dim(data)[2]
  Kn <- array(0, c(n))
  u.upper <- array(0, c(n, d))
  u.lower <- array(0, c(n, d))
  (1/n) * (3/4)^d * sum(sapply(1:nrow(data), 
    function(k) {
      u.upper <- (x + l - data[k,]) / h
      if (any(u.upper < -1)) 0
      else {
        u.lower <- (x - data[k,]) / h
        if (any(u.lower > 1)) 0
        else {      
          u.upper[u.upper > 1] <- 1
          u.lower[u.lower < -1] <- -1            
          prod(u.upper - 1/3 * u.upper^3 - u.lower + 1/3 * u.lower^3)    
        }
      }
    }
  )) 
  
}


d <- c(2, 4, 4, 4, 4, 10, 12, 12, 15, 16, 20, 50)
d <- as.matrix(d)

d.ecdf <- ecdf(d)
plot(d.ecdf, xlim=c(0,max(d)))


h1 <- c(0.001)
h2 <- c(0.0001, 0.0001)
h3 <- c(3)
h4 <- c(5)

x <- c(4.5)

len = 100

a <- seq(0, max(d), length.out=len)
or <- rep(0, len)
V1 <- array(0, c(len))
V2 <- array(0, c(len))
V3 <- array(0, c(len))
V4 <- array(0, c(len))

for (i in 1:len) {
  #V1[i] <- H(d, c(a[i]), h1)
  #lines(a, V1, col="red", type="l")
  
  #V2[i] <- H(d, c(a[i]), h2)
  #lines(a,V2, col="green")
  
  #V3[i] <- H(d, c(a[i]), h3)
  #lines(a, V3, col="black")
  
  V4[i] <- H2(d, c(or[i]), c(a[i]), h2)    
}
lines(a, V4, col="blue", type="l")



